##### This script runs the linear mixed-effects models reported in the Table OA2.

# The following set of models do not control for compositional effects
# (because the individual-level variables are centered within clusters CWC)

# Model A: two levels with RE at the COUNTRY-ROUND level
Model_A_CWC <- lmer(red_sup ~ employed_CWC + unemployed_CWC + hinc_2_cop_CWC + hinc_3_dif_CWC + hinc_4_vdif_CWC + female_CWC + age_CWC + isced2_CWC + isced3_CWC + isced4_CWC + isced56_CWC + stfeco_CWC + lrscale_CWC + 
                  LFPR_y0 + SOCEXP_y0 + GDP_y0 + gini_disp_y0 + liberal + mediter + (1 | cntry_round), data = ess_4_ld_models, weights = pspwght)
summary(Model_A_CWC)
check_collinearity(Model_A_CWC)


# Model B: two levels with RE at the COUNTRY level
Model_B_CWC <- lmer(red_sup ~ employed_CWC + unemployed_CWC + hinc_2_cop_CWC + hinc_3_dif_CWC + hinc_4_vdif_CWC + female_CWC + age_CWC + isced2_CWC + isced3_CWC + isced4_CWC + isced56_CWC + stfeco_CWC + lrscale_CWC + 
                  LFPR_y0 + SOCEXP_y0 + GDP_y0 + gini_disp_y0 + liberal + mediter + (1 | cntry), data = ess_4_ld_models, weights = pspwght)
summary(Model_B_CWC)
check_collinearity(Model_B_CWC)


# Model C: three levels with RE at ROUND level (3) and then COUNTRY-ROUND level (2)
Model_C_CWC <- lmer(red_sup ~ employed_CWC + unemployed_CWC + hinc_2_cop_CWC + hinc_3_dif_CWC + hinc_4_vdif_CWC + female_CWC + age_CWC + isced2_CWC + isced3_CWC + isced4_CWC + isced56_CWC + stfeco_CWC + lrscale_CWC + 
                  LFPR_y0 + SOCEXP_y0 + GDP_y0 + gini_disp_y0 + liberal + mediter + (1 | cntry_round) + (1 | essround), data = ess_4_ld_models, weights = pspwght)
summary(Model_C_CWC)
check_collinearity(Model_C_CWC)


# Model D: three levels with RE at COUNTRY level (3) and then COUNTRY-ROUND level (2)
Model_D_CWC <- lmer(red_sup ~ employed_CWC + unemployed_CWC + hinc_2_cop_CWC + hinc_3_dif_CWC + hinc_4_vdif_CWC + female_CWC + age_CWC + isced2_CWC + isced3_CWC + isced4_CWC + isced56_CWC + stfeco_CWC + lrscale_CWC + 
                  LFPR_y0 + SOCEXP_y0 + GDP_y0 + gini_disp_y0 + liberal + mediter + (1 | cntry_round) + (1 | cntry), data = ess_4_ld_models, weights = pspwght)
summary(Model_D_CWC)
check_collinearity(Model_D_CWC)


# Model E: cross-classified model
Model_E_CWC <- lmer(red_sup ~ employed_CWC + unemployed_CWC + hinc_2_cop_CWC + hinc_3_dif_CWC + hinc_4_vdif_CWC + female_CWC + age_CWC + isced2_CWC + isced3_CWC + isced4_CWC + isced56_CWC + stfeco_CWC + lrscale_CWC + 
                  LFPR_y0 + SOCEXP_y0 + GDP_y0 + gini_disp_y0 + liberal + mediter + (1 | essround) + (1 | cntry), data = ess_4_ld_models, weights = pspwght)
summary(Model_E_CWC)
check_collinearity(Model_E_CWC)


# Model F: full model
Model_F_CWC <- lmer(red_sup ~ employed_CWC + unemployed_CWC + hinc_2_cop_CWC + hinc_3_dif_CWC + hinc_4_vdif_CWC + female_CWC + age_CWC + isced2_CWC + isced3_CWC + isced4_CWC + isced56_CWC + stfeco_CWC + lrscale_CWC + 
                  LFPR_y0 + SOCEXP_y0 + GDP_y0 + gini_disp_y0 + liberal + mediter + (1 | cntry_round) + (1 | essround) + (1 | cntry), data = ess_4_ld_models, weights = pspwght)
summary(Model_F_CWC)
check_collinearity(Model_F_CWC)


# Model G: two levels with RE at the YEAR level
Model_G_CWC <- lmer(red_sup ~ employed_CWC + unemployed_CWC + hinc_2_cop_CWC + hinc_3_dif_CWC + hinc_4_vdif_CWC + female_CWC + age_CWC + isced2_CWC + isced3_CWC + isced4_CWC + isced56_CWC + stfeco_CWC + lrscale_CWC + 
                  LFPR_y0 + SOCEXP_y0 + GDP_y0 + gini_disp_y0 + liberal + mediter + (1 | essround), data = ess_4_ld_models, weights = pspwght)
summary(Model_G_CWC)
check_collinearity(Model_G_CWC)


# Plotting the model results
tab_model(Model_A_CWC, Model_B_CWC, Model_C_CWC, Model_D_CWC, Model_E_CWC, Model_F_CWC, Model_G_CWC, digits = 3, digits.p = 4, digits.re = 4)



#################################################################################
# Storing and exporting results of the seven models into a Table OA2 of the paper #
#################################################################################

# Creating an empty data.frame into which the results will be stored
Models_CWC_individual <- data.frame(c("Intercept", "", "employed_CWC", "", "unemployed_CWC", "", 
                               "hinc_2_cop_CWC", "", "hinc_3_dif_CWC", "", "hinc_4_vdif_CWC", "", 
                               "female_CWC", "", "age_CWC", "", "isced2_CWC", "", "isced3_CWC", "", 
                               "isced4_CWC", "", "isced56_CWC", "", "stfeco_CWC", "", "lrscale_CWC", "", 
                               "LFPR_y0", "", "SOCEXP_y0", "", "GDP_y0", "", "gini_disp_y0", "", 
                               "liberal", "", "mediter", "", "Var_cntry", "Var_essround", 
                               "Var_cntry_round", "Var_individual", "ICC_cntry", "ICC_essround", "ICC_cntry_round"), matrix(data = rep(NA, 47*7), nrow = 47))
colnames(Models_CWC_individual) <- c("Parameter", "Model A", "Model B", "Model C", "Model D", "Model E", "Model F", "Model G")

# Using the programmed function assign_fe_results to create Table OA2 of the article
Models_CWC_individual <- assign_fe_results(model_name = Model_A_CWC, assignment_object = Models_CWC_individual, column_number = 2)
Models_CWC_individual <- assign_fe_results(model_name = Model_B_CWC, assignment_object = Models_CWC_individual, column_number = 3)
Models_CWC_individual <- assign_fe_results(model_name = Model_C_CWC, assignment_object = Models_CWC_individual, column_number = 4)
Models_CWC_individual <- assign_fe_results(model_name = Model_D_CWC, assignment_object = Models_CWC_individual, column_number = 5)
Models_CWC_individual <- assign_fe_results(model_name = Model_E_CWC, assignment_object = Models_CWC_individual, column_number = 6)
Models_CWC_individual <- assign_fe_results(model_name = Model_F_CWC, assignment_object = Models_CWC_individual, column_number = 7)
Models_CWC_individual <- assign_fe_results(model_name = Model_G_CWC, assignment_object = Models_CWC_individual, column_number = 8)
View(Models_CWC_individual)

### Manual assignment of variance components and ICCs
# Model A: RE at cntry_round
Models_CWC_individual[Models_CWC_individual$Parameter == "Var_cntry_round", "Model A"] <- as.data.frame(VarCorr(Model_A_CWC))[, "vcov"][1]
Models_CWC_individual[Models_CWC_individual$Parameter == "Var_individual", "Model A"] <- as.data.frame(VarCorr(Model_A_CWC))[, "vcov"][2]
icc(Model_A_CWC, by_group = TRUE)
Models_CWC_individual[Models_CWC_individual$Parameter == "ICC_cntry_round", "Model A"] <- (as.data.frame(VarCorr(Model_A_CWC))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_A_CWC))[, "vcov"]))*100


# Model B: RE at cntry
Models_CWC_individual[Models_CWC_individual$Parameter == "Var_cntry", "Model B"] <- as.data.frame(VarCorr(Model_B_CWC))[, "vcov"][1]
Models_CWC_individual[Models_CWC_individual$Parameter == "Var_individual", "Model B"] <- as.data.frame(VarCorr(Model_B_CWC))[, "vcov"][2]
icc(Model_B_CWC, by_group = TRUE)
Models_CWC_individual[Models_CWC_individual$Parameter == "ICC_cntry", "Model B"] <- (as.data.frame(VarCorr(Model_B_CWC))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_B_CWC))[, "vcov"]))*100

# Model C: RE at essround and cntry_round
Models_CWC_individual[Models_CWC_individual$Parameter == "Var_essround", "Model C"] <- as.data.frame(VarCorr(Model_C_CWC))[, "vcov"][2]
Models_CWC_individual[Models_CWC_individual$Parameter == "Var_cntry_round", "Model C"] <- as.data.frame(VarCorr(Model_C_CWC))[, "vcov"][1]
Models_CWC_individual[Models_CWC_individual$Parameter == "Var_individual", "Model C"] <- as.data.frame(VarCorr(Model_C_CWC))[, "vcov"][3]
icc(Model_C_CWC, by_group = TRUE)
Models_CWC_individual[Models_CWC_individual$Parameter == "ICC_essround", "Model C"] <- (as.data.frame(VarCorr(Model_C_CWC))[, "vcov"][2]/sum(as.data.frame(VarCorr(Model_C_CWC))[, "vcov"]))*100
Models_CWC_individual[Models_CWC_individual$Parameter == "ICC_cntry_round", "Model C"] <- (as.data.frame(VarCorr(Model_C_CWC))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_C_CWC))[, "vcov"]))*100

# Model D: three levels with RE at COUNTRY level (3) and then COUNTRY-ROUND level (2)
Models_CWC_individual[Models_CWC_individual$Parameter == "Var_cntry", "Model D"] <- as.data.frame(VarCorr(Model_D_CWC))[, "vcov"][2]
Models_CWC_individual[Models_CWC_individual$Parameter == "Var_cntry_round", "Model D"] <- as.data.frame(VarCorr(Model_D_CWC))[, "vcov"][1]
Models_CWC_individual[Models_CWC_individual$Parameter == "Var_individual", "Model D"] <- as.data.frame(VarCorr(Model_D_CWC))[, "vcov"][3]
icc(Model_D_CWC, by_group = TRUE)
Models_CWC_individual[Models_CWC_individual$Parameter == "ICC_cntry", "Model D"] <- (as.data.frame(VarCorr(Model_D_CWC))[, "vcov"][2]/sum(as.data.frame(VarCorr(Model_D_CWC))[, "vcov"]))*100
Models_CWC_individual[Models_CWC_individual$Parameter == "ICC_cntry_round", "Model D"] <- (as.data.frame(VarCorr(Model_D_CWC))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_D_CWC))[, "vcov"]))*100

# Model E: RE at country and ess_round
Models_CWC_individual[Models_CWC_individual$Parameter == "Var_cntry", "Model E"] <- as.data.frame(VarCorr(Model_E_CWC))[, "vcov"][1]
Models_CWC_individual[Models_CWC_individual$Parameter == "Var_essround", "Model E"] <- as.data.frame(VarCorr(Model_E_CWC))[, "vcov"][2]
Models_CWC_individual[Models_CWC_individual$Parameter == "Var_individual", "Model E"] <- as.data.frame(VarCorr(Model_E_CWC))[, "vcov"][3]
icc(Model_E_CWC, by_group = TRUE)
Models_CWC_individual[Models_CWC_individual$Parameter == "ICC_cntry", "Model E"] <- (as.data.frame(VarCorr(Model_E_CWC))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_E_CWC))[, "vcov"]))*100
Models_CWC_individual[Models_CWC_individual$Parameter == "ICC_essround", "Model E"] <- (as.data.frame(VarCorr(Model_E_CWC))[, "vcov"][2]/sum(as.data.frame(VarCorr(Model_E_CWC))[, "vcov"]))*100

# Model F: RE at all three effects
Models_CWC_individual[Models_CWC_individual$Parameter == "Var_cntry", "Model F"] <- as.data.frame(VarCorr(Model_F_CWC))[, "vcov"][2]
Models_CWC_individual[Models_CWC_individual$Parameter == "Var_essround", "Model F"] <- as.data.frame(VarCorr(Model_F_CWC))[, "vcov"][3]
Models_CWC_individual[Models_CWC_individual$Parameter == "Var_cntry_round", "Model F"] <- as.data.frame(VarCorr(Model_F_CWC))[, "vcov"][1]
Models_CWC_individual[Models_CWC_individual$Parameter == "Var_individual", "Model F"] <- as.data.frame(VarCorr(Model_F_CWC))[, "vcov"][4]
icc(Model_F_CWC, by_group = TRUE)
Models_CWC_individual[Models_CWC_individual$Parameter == "ICC_cntry", "Model F"] <- (as.data.frame(VarCorr(Model_F_CWC))[, "vcov"][2]/sum(as.data.frame(VarCorr(Model_F_CWC))[, "vcov"]))*100
Models_CWC_individual[Models_CWC_individual$Parameter == "ICC_essround", "Model F"] <- (as.data.frame(VarCorr(Model_F_CWC))[, "vcov"][3]/sum(as.data.frame(VarCorr(Model_F_CWC))[, "vcov"]))*100
Models_CWC_individual[Models_CWC_individual$Parameter == "ICC_cntry_round", "Model F"] <- (as.data.frame(VarCorr(Model_F_CWC))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_F_CWC))[, "vcov"]))*100

# Model G: RE at essround
Models_CWC_individual[Models_CWC_individual$Parameter == "Var_essround", "Model G"] <- as.data.frame(VarCorr(Model_G_CWC))[, "vcov"][1]
Models_CWC_individual[Models_CWC_individual$Parameter == "Var_individual", "Model G"] <- as.data.frame(VarCorr(Model_G_CWC))[, "vcov"][2]
icc(Model_G_CWC, by_group = TRUE)
Models_CWC_individual[Models_CWC_individual$Parameter == "ICC_essround", "Model G"] <- (as.data.frame(VarCorr(Model_G_CWC))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_G_CWC))[, "vcov"]))*100

# Exporting Table OA2 to a csv file
write.csv(x = Models_CWC_individual, file = "table_OA2_Models_CWC_individual.csv")



### Collecting the coefficient estimates for the LFPR contextual variable and storing them to a data.frame
LFPR_CWC_individual <- data.frame(c("Model A", "Model B", "Model C", "Model D", "Model E", "Model F", "Model G"),
                                  matrix(data = rep(NA, 7*3), nrow = 7))
colnames(LFPR_CWC_individual) <- c("Model type", "Estimate", "Std. Error", "df")

LFPR_CWC_individual[LFPR_CWC_individual$`Model type` == "Model A", 2:4] <- coef(summary(Model_A_CWC))[15, c("Estimate", "Std. Error", "df")]
LFPR_CWC_individual[LFPR_CWC_individual$`Model type` == "Model B", 2:4] <- coef(summary(Model_B_CWC))[15, c("Estimate", "Std. Error", "df")]
LFPR_CWC_individual[LFPR_CWC_individual$`Model type` == "Model C", 2:4] <- coef(summary(Model_C_CWC))[15, c("Estimate", "Std. Error", "df")]
LFPR_CWC_individual[LFPR_CWC_individual$`Model type` == "Model D", 2:4] <- coef(summary(Model_D_CWC))[15, c("Estimate", "Std. Error", "df")]
LFPR_CWC_individual[LFPR_CWC_individual$`Model type` == "Model E", 2:4] <- coef(summary(Model_E_CWC))[15, c("Estimate", "Std. Error", "df")]
LFPR_CWC_individual[LFPR_CWC_individual$`Model type` == "Model F", 2:4] <- coef(summary(Model_F_CWC))[15, c("Estimate", "Std. Error", "df")]
LFPR_CWC_individual[LFPR_CWC_individual$`Model type` == "Model G", 2:4] <- coef(summary(Model_G_CWC))[15, c("Estimate", "Std. Error", "df")]

# Computing the 95% confidence interval (based on student?s t-distribution)
LFPR_CWC_individual$lower_bound <- LFPR_CWC_individual$Estimate - LFPR_CWC_individual$`Std. Error`*qt(p = 0.975, df = LFPR_CWC_individual$df)
LFPR_CWC_individual$upper_bound <- LFPR_CWC_individual$Estimate + LFPR_CWC_individual$`Std. Error`*qt(p = 0.975, df = LFPR_CWC_individual$df)



# Deleting model results from the workspace (as they are no longer needed)
rm(Model_A_CWC, Model_B_CWC, Model_C_CWC, Model_D_CWC, Model_E_CWC, Model_F_CWC, Model_G_CWC)
